{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Viewing and manipulating FITS images\n",
    "\n",
    "## Authors\n",
    "Lia Corrales, Kris Stern, Stephanie T. Douglas, Kelle Cruz\n",
    "\n",
    "## Learning Goals\n",
    "* Open FITS files and load image data\n",
    "2. Make a 2D histogram with image data\n",
    "3. Stack several images into a single image\n",
    "4. Write image data to a FITS file\n",
    "\n",
    "## Keywords\n",
    "FITS, file input/output, image manipulation, numpy, matplotlib, histogram, colorbar\n",
    "\n",
    "## Summary\n",
    "\n",
    "This tutorial demonstrates the use of `astropy.utils.data` to download a data file, then uses `astropy.io.fits` to open the file, and lastly uses `matplotlib` to view the image with different color scales and stretches and to make histograms. In this tutorial we've also included a demonstration of simple image stacking."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Set up matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "from astropy.io import fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download the example FITS files for this tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.utils.data import download_file\n",
    "image_file = download_file('http://data.astropy.org/tutorials/FITS-images/HorseHead.fits', cache=True )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Opening FITS files and loading the image data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's open the FITS file to find out what it contains."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hdu_list = fits.open(image_file)\n",
    "hdu_list.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generally, the image information is located in the <code>PRIMARY</code> block. The blocks are numbered and can be accessed by indexing <code>hdu_list</code>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image_data = hdu_list[0].data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our data is now stored as a 2D numpy array.  But how do we know the dimensions of the image?  We can look at the `shape` of the array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(type(image_data))\n",
    "print(image_data.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great! At this point, we can close the FITS file because we've stored everything we wanted to a variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hdu_list.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SHORTCUT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you don't need to examine the FITS header, you can call `fits.getdata` to bypass the previous steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image_data = fits.getdata(image_file)\n",
    "print(type(image_data))\n",
    "print(image_data.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Viewing the image data and getting basic statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(image_data, cmap='gray')\n",
    "plt.colorbar()\n",
    "\n",
    "# To see more color maps\n",
    "# http://wiki.scipy.org/Cookbook/Matplotlib/Show_colormaps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's get some basic statistics about our image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Min:', np.min(image_data))\n",
    "print('Max:', np.max(image_data))\n",
    "print('Mean:', np.mean(image_data))\n",
    "print('Stdev:', np.std(image_data))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting a histogram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make a histogram with `matplotlib.pyplot.hist()`, we'll need to cast the data from a 2D array to something one dimensional."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, let's use the `ndarray.flatten()` to return a 1D numpy array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(type(image_data.flatten()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "histogram = plt.hist(image_data.flatten(), bins='auto')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Displaying the image with a logarithmic scale"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What if we want to use a logarithmic color scale? To do so, we'll need to load the `LogNorm` object from `matplotlib`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.colors import LogNorm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(image_data, cmap='gray', norm=LogNorm())\n",
    "\n",
    "# I chose the tick marks based on the histogram above\n",
    "cbar = plt.colorbar(ticks=[5.e3,1.e4,2.e4])\n",
    "cbar.ax.set_yticklabels(['5,000','10,000','20,000'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic image math: image stacking"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also perform math with the image data like any other numpy array.  In this particular example, we'll stack several images of M13 taken with a ~10'' telescope."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start by opening a series of FITS files and store the data in a list, which we've named `image_concat`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base_url = 'http://data.astropy.org/tutorials/FITS-images/M13_blue_{0:04d}.fits'\n",
    "\n",
    "image_list = [download_file(base_url.format(n), cache=True) \n",
    "              for n in range(1, 5+1)]\n",
    "image_concat = [fits.getdata(image) for image in image_list]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll stack the images by summing the concatenated list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The long way\n",
    "final_image = np.zeros(shape=image_concat[0].shape)\n",
    "\n",
    "for image in image_concat:\n",
    "    final_image += image\n",
    "\n",
    "# The short way\n",
    "# final_image = np.sum(image_concat, axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to show the image, but need to decide on the best stretch. To do so let's plot a histogram of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image_hist = plt.hist(final_image.flatten(), bins='auto')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use the keywords `vmin` and `vmax` to set limits on the color scaling for `imshow`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(final_image, cmap='gray', vmin=2E3, vmax=3E3)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Writing image data to a FITS file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can do this with the `writeto()` method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Warning: you'll receive an error if the file you are trying to write already exists.  That's why we've set `overwrite=True`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outfile = 'stacked_M13_blue.fits'\n",
    "\n",
    "hdu = fits.PrimaryHDU(final_image)\n",
    "hdu.writeto(outfile, overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determine the mean, median, and standard deviation of a part of the stacked M13 image where there is *not* light from M13.  Use those statistics with a sum over the part of the image that includes M13 to estimate the total light in this image from M13."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show the image of the Horsehead Nebula, but in units of *surface brightness* (magnitudes per square arcsecond).\n",
    "(Hint: the *physical* size of the image is 15x15 arcminutes.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now write out the image you just created, preserving the header the original image had, but add a keyword 'UNITS' with the value 'mag per sq arcsec'.\n",
    "(Hint: it may be helpful to read the [astropy.io.fits](http://docs.astropy.org/en/stable/io/fits/index.html) documentation if you're not sure how to include both the header and the data.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "astropy-tutorials": {
   "author": "Lia R. Corrales <lia@astro.columbia.edu>",
   "date": "January 2014",
   "description": "Demonstrates use of astropy.utils.data to download data; astropy.io.fits to open the file; matplotlib to view the image with different color scales and stretches and to make histrograms. Also includes method for simple image stacking.",
   "link_name": "Viewing and manipulating FITS images",
   "name": "",
   "published": true
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
